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We solve the Gross-Pitaevskii equation to model the behaviour of a weakly-interacting 
Bose condensate. Solutions are presented for the eigenstates in one- and two-dimensional 
harmonic traps. We include the effect of gravity and coupling to a second condensate 
to simulate an output coupler for Bose condensed atoms, and find that the output pulse 
shape is in qualitative agreement with experiment. We also model flow under gravity 
through a constriction, demonstrating excitation and separation of eigenmodes in the 
condensate. 



I. INTRODUCTION 

The first experimental observations of Bose-Einstein condensation in magnetically-trapped alkali atoms 
0-^1 excited considerable interest in the properties of weakly-interacting Bose fluids. Subsequent exper- 
iments have studied collective excitations and sound propagation in the condensate ^,0], leading to 
hopes of a deeper understanding of the nature of superfluidity B . An important recent experiment has 
shown macroscopic interference between two expanding condensates, providing compelling evidence for 
coherence in the ground state of the trap M. Thus, the demonstration of an output coupler for Bose 



condensed atoms |10 may be thought of as realizing a simple pulsed atom laser, though controversy still 
surrounds this issue p"l| ]. 

Theoretical work on weakly-interacting Bose condensates has focused on solutions of the Gross-Pitaevskii 
or time-dependent nonlinear Schrodinger equation (NLSE) p"^-|i~5|] . This equation describes the evolution 
of a self-consistent wavefunction representing a dilute near-zero temperature condensate, and follows from 
mean- field theory where interactions between atoms are modelled by an s-wave scattering length 
The Gross-Pitaevskii equation has been shown to provide an accurate description of the shape |ll| and 
dynamics p|]5|]l7| of weakly interacting condensates, with the principal effect of finite temperature arising 
as condensate depletion Jig) . 

The time-dependent NLSE may be solved numerically by a variety of methods Jll| . Previous studies have 
used semi- implicit Crank-Nicholson in one dimension jl3| |, or alternating-direction implicit-difference in 
two dimensions [jl5). In this work, we solve the NLSE using Split-Step Fourier Jll],|2(J and finite difference 
methods. A combination of these techniques can be used to model the dynamics of a condensate under 
a wide range of situations, for example, the creation of vortices by a moving object 
Here we consider two examples: simulations of the output coupling of Bose condensed atoms from a trap; 
and the flow of a condensate through a constriction. The output coupling of Bose condensed atoms, 
by driving an RF transition to a untrapped state, has been demonstrated experimentally by Mewes 
et al. ]i0[. Previous theoretical studies have considered only ID motion without gravity p2|, or have 



modelled the process by neglecting kinetic energy [£3[ or interaction terms |24j in the Gross-Pitaevskii 
equation. We solve coupled Gross-Pitaevskii equations in 2D, with the inclusion of gravity, to model the 
evolution of both the output pulse and the remaining trapped condensate. The simulations reproduce the 
crescent-shaped atomic pulses observed experimentally Jl(],^5| , and we demonstrate that the recoil of the 
output pulse excites a harmonic oscillation of the parent condensate. Secondly, we consider techniques 
for condensate manipulation. To maintain the shape of a condensate pulse, it could be coupled into 
a waveguide such as a collimated far-off resonance laser beam. However, a more interesting situation 
arises in the flow through a constriction. For this case, we predict the excitation and spatial filtering of 
eigenmodes. 

The paper is arranged as follows: Section 2 contrasts the Fourier and finite-difference methods, and 
discusses ID and 2D solutions of the Gross-Pitaevskii equation for the lowest energy levels. Simulations 
of the output coupling and flow through a constriction are presented in Sections 3 and 4 respectively. 
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II. NUMERICAL SOLUTION OF THE GROSS-PITAEVSKII EQUATION 



The property of phase coherence inherent in Bose condensation implies that it may be represented by an 
order parameter, which is simply the macroscopic wavefunction of the condensate. At zero temperature 
this wavefunction, ^(r,t), is described by the Gross-Pitaevskii equation. Including a trap potential, 
V(r,t), this has the form |||: 

ih^(r,t)= (-^V 2 + V(r,t)+NU \*(r,t)\ 2 ^<f(r,t), (1) 

where m is the atomic mass, N is the number of atoms, and Uq describes the interaction between atoms 
in the condensate. For T w interactions are modelled by the s-wave scattering length, a, such that 
Uq = Airh 2 a/m. 

First, consider a ID condensate, where the trapping potential is described by V(x) = muj 2 x 2 /2. For 
convenience, equation (Q) is scaled in terms of harmonic oscillator units (h.o.u.) of length and time, where 
£ = (h/2muj)~ 1 / 2 x and t = cot. With the change in length scale the normalisation of the wavefunction is 
also modified, and thus t/'(£j t ) = (fr/2rnuj) 1 / 4:1 § l (x 1 t). Taking the unit of energy to be hu>, equation ([!]) 
becomes 

i§ ? mr)=(-^ + \e + C\mr)\ 2 ^mr), (2) 

where C = AitNa^Ti/moj) 1 / 2 . In the next two subsections, we discuss the solution of this NLSE using 
(A) the Split-Step Fourier method, and (B) finite differences. 



A. Fourier method 



Consider a Schrodinger equation of the form: 

ih^-Mr, t) = [T+ V(r, t)ty(r, t) = H{r, i)^(r, t), (3) 
at 

with the solution, ip(r,0), at t = 0. We can propagate the solution through a short time interval to 
second-order acccuracy using the unitary evolution operator. Expansion of the time-evolution operator 
gives §|: 



-iHt/H _ e -iV(0)t/2h e ~iTt/h e -iV(0)t/2H + Q(t 2 ). (4) 



For a slowly- varying V(t), the order of accuracy of this 'half-step' expansion increases to 0(t ). The 



operator, exp(—iTt/h) can be evaluated numerically by a variety of methods |19|. The Split-Step Fourier 
method relies on the use of discrete Fourier transforms computed with a Fast Fourier Transform (FFT) 
algorithm, allowing high computational efficiency. Moreover, the Split-Step Fourier method is readily 
extended to solution of the NLSE in an arbitrary number of dimensions by use of iV-dimensional FFT 
routines p7| . 

Solutions of (|^) may be evaluated numerically by beginning with an analytic ground or excited state 
solution in the absence of the nonlinear term. The solution is then propagated through real time using 
the Split-Step Fourier method, while at each time step, the value of the nonlinear constant C is increased 
adiabatically ('ramped') until reaching the desired value of C . In practice, this consists of inserting a time- 
dependent pre-factor, P(t) — [1 — cos(7rr/r r )]/2, in front of the nonlinear term fl5|| , where the ramp time 
r r is chosen to be long compared to the oscillation time of the ground state of the trap. Alternatively, the 
ground state can be readily determined through propagation in imaginary time (t — > —it) pcfl . However, 
we find that imaginary time FFT is inappropriate for the excited states, as these are exponentially 
suppressed with respect to lower energy modes. 

One-dimensional wavefunctions for the first and second excited states are presented in Figs. 1 and 2, and 
show interesting behaviour as the repulsive interactions increase in strength. In both cases, the edges of 
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the wavefunction spread out in a similar fashion to the ground state J13| . However, for the first excited 
state the width of the central node barely changes as C increases; while in the second excited state, the 
two nodes of the symmetric condensate asymptotically approach a fixed separation, A£ ~ 1, for large 
C. Later, we will demonstrate how such modes could be excited by flow of a condensate through a 
constriction (Section 4). 
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FIG. 1. Condensate density \i/j\ 2 plotted against position (in h.o.u.) for the first excited state in ID, with the 
nonlinear coefficient C = 0, 30, 60, 90, 120. As C increases, the edges of the symmetric condensate spread out 
while the central node remains almost unchanged. 
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FIG. 2. Condensate density \ifj\ 2 plotted against position (in h.o.u.) for the second excited state in ID, for 
C — 0, 30, 60, 90, 120, 150. As the interactions increase the nodes of the symmetric condensate converge, and tend 
towards a fixed separation at high C. 
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B. Finite differences 



The procedure described above is difficult to implement for C > 1 50, requiring long computation times — 
a problem compounded in 2D. An alternative approach is to evaluate the solutions of the NLSE and 
corresponding chemical potentials using finite differences. By replacing the differentials in @ by finite 
differences between all adjacent grid points, and including the condition that the wavefunction is nor- 
malised, we obtain a set of simultaneous equations for values of the wavefunction at each grid point. 
The simultaneous nonlinear equations can then be solved by iteration, e.g. using the Newton-Raphson 
method. Table I displays chemical potentials for the ground and first two excited states of the trap, in 
units of hui. By neglecting the kinetic term in the time- independent NLSE — an approximation valid for 
large interaction strengths — one obtains a simple analytic solution for the ground state |]l2"f , for which the 
chemical potential is [1\d = (3C/8) 2 / 3 . This is commonly known as the Thomas-Fermi approximation. 
At large C, the data for the ground state does indeed approach the Thomas-Fermi limit. In addition, 
the ground and excited energies converge so that they approach a fixed, equal separation at high C . 
For a 2D isotropic trap, V(x,z) — mcu 2 (x 2 + z 2 )/2, we can reduce the problem to ID by converting to 
cylindrical polar coordinates. Then the eigenvalue equation has the form: 

- \h - 7W + i" 2 - " + CWp - "• t)|2 ) *<" » T) = (5) 

where the ground state (i.e. no (^-dependence) is to be found. Ground-state chemical potentials for the 2D 
circular trap are presented in Table I. An important point to note is that the results approach the Thomas- 
Fermi limit, [i2D — (C/2tt) 1 ^ 2 , more slowly than in ID. The parabolic Thomas-Fermi wavefunction fails 
to match the exact solution at the edges of the condensate. In ID the wavefunctions mismatch at only 
two points, while in 2D the mismatched region is a ring. Thus, as the dimensionality of the condensate 
increases, the Thomas-Fermi solution attains less accuracy, and one must be cautious when using this 
approximation in 3D. 

This example also illustrates the ease with which the finite-difference method can find solutions of a 
NLSE containing large nonlinearities, when compared to the Split-Step Fourier method. However, our 
experience, in agreement with detailed studies fll9| , indicates that the Fourier method is more efficient 
when propagating the solution subject to a time-dependent potential. This is of particular significance 
when considering 2D simulations, where it is important that the computational times remain manageable. 
So, in the 2D simulations contained in this paper, we first find the ground-state solution of the NLSE for 
a trapped condensate using the finite difference method, and then propagate the wavefunction under a 
time-dependent potential using the Split-Step Fourier method. 



TABLE I. Chemical potentials (in units of hui) for the three lowest energy states in ID, and for the ground 
state of an isotropic 2D condensate, calculated using finite differences. Ground state values are compared to the 
corresponding Thomas-Fermi values: /iir> = (3C/8) 2//3 and fi2D = (C /2-k) 1 ^ 2 . 
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III. THE OUTPUT COUPLER 



A. The coupled Gross-Pitaevskii equations 



The previous section discussed how a combination of Fourier and finite difference methods are particularly 
efficient at finding and propagating solutions of the Gross-Pitaevskii equation. The same technique can be 
applied to mixtures of condensates. In particular, we solve coupled equations to describe the experimental 
realization of an output coupler for Bose condensed atoms, where a sequence of 'daughter' condensate 
pulses is produced by excitation of a 'parent' condensate 119]. p or 23 Na (or 87 Rb Jl[) condensates, the 
F = 1 ground state is composed of the trapped state Mp — —1, the magnetically insensitive level 
Mp = 0, and the anti-trapped state Mp = 1. Condensed atoms in each state may be represented 
by a Bose wavefunction ^(r,*), where i 6 { — 1,0, +1}. The Hamiltonian of the system can then be 
constructed, to include gravitational and magnetic trapping potentials, hyperfine interactions, and the 
RF coupling term — fi ■ B(r,t). Our approach involves applying the variational principle to obtain the 
following equations describing the dynamics of the parent (^-i) and daughter condensates (\I/o,^i) : 



d 

ih—V k (r,t) 



2m * — ' 



tffctr.tj + ftj^fijytf^r.t), 



(6) 



where Uk(r) represents the gravitational and trap potentials in the substate k, and Q k j is the Rabi 
frequency describing coupling between sublevels. Atom-atom scattering appears as the matrix element 
(Wjlyl^j). Assuming all scattering lengths to be equal, V may be replaced by the pseudopotential 
U 5(ri -rj). 

In the experiment, the RF pulse was of short duration and small area JTcj] , leading to the following 
simplifications. First, due to the small area (| / Sl+ifldt\ <C 1), coupling to the Mp = — 1 state is small 
and may be safely ignored, and also back-coupling, i.e. Rabi oscillations, are negligible. An example 
of strong coupling, where Rabi oscillations are important, has been considered by Ballagh et al. p^] in 
a ID weightless model. Secondly, the short pulse duration, T results in broadband excitation, so that 
we may assume resonant coupling across the trap and neglect the spatial dependence of the coupling 
term, fi(r,i'). In addition, we may assume that the parent condensate is 'frozen' during the interaction 
(uiT «C 1), in which case the daughter condensate is created suddenly as a shadow of the parent, i.e.: 



l*o(r,T)| 



dt / n(r,* / )*-i(r,0 «/|*-i(r,0)| 



(7) 



and, 



|*U(r,T)| a «(i-/)|iiu(r,o)| 5 



(8) 



where /, the fraction of atoms coupled out of the trap, is equal to the square of the pulse area. 
Following these simplifications, we obtain a pair of coupled time-dependent equations in 2D: 

idrif-i = (-v 2 + \(e + v 2 ) + c\^-i\ 2 + c\M 2 ) V'-i, 



(9) 



id T ^o = (-V 2 + Gr, + C\i> \ 2 + C|V-i| 2 ) Vo, (10) 

where £ = (h/2muj)- 1/2 x, T] = (fi/2mw)^/ 2 z, and ipi(£,ri,T) = (li/2mw)" 1/2 $ t (j;, z, t). Thus, C = 
8>TtNa, where N is the linear density of atoms in the y-direction, and gravity is parametrized by G = 
(m/2froj 3 ) 1 / 2 g. Note that gravity is absent from Equation @, as this term only gives rise to displacement 
of the equilibrium position of the trapped condensate. With initial conditions given by Equations (0) 
and (^) , equations (Q) and ( |To| ) are solved to model the dynamical behaviour of the parent and daughter 
condensates for r > uT. The results of these simulations are now discussed. 
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B. Results 



To simulate the MIT experiment we consider the parameters C — 200 and G — 3.0, corresponding to a 
trap frequency of ui ~ 2tt x 200 Hz and N « 2 x 10 4 23 Na atoms per h.o.u. respectively. Density plots 
of the output pulse (for / = 0.2) are presented in Fig. 3. After creation, the output pulse is subject to 
gravity, a repulsive internal force, and a strong repulsion from the residual parent condensate. We see 
in Fig. 3 that even at early times, a combination of these effects dramatically changes the shape of the 
output. The curvature of the daughter condensate is found to be less pronounced for weaker interaction 
strengths. At later times the output pulse expands freely, resulting in a 'crescent-shaped' density profile, 
characteristic of the MIT experiment This suggests that the zero-temperature approximation, 

inherent in the coupled Gross-Pitaevskii equations, is appropriate for modelling the dynamics of the 
output. The curvature of the output pulse was not observed in previous theoretical studies [ p4| , due to 
a different trap geometry and subsequent neglect of the interaction term in the axial direction. 




FIG. 3. Contour plot of \ip(£, Vy r )l f° r an output coupled Bose condensate falling under gravity, with C = 200, 
/ = 0.20 and G = 3.0, at times r = (offset f -> $-30), r = 1.2 (offset £ -» £-10), and r = 2.4 (offset £ -> £ + 20). 
Nine equally-spaced contours are drawn for each time-frame, up to peak densities of |i/j(0,0,0)| 2 = 5.67 x 10 -3 , 
|^(0, -7.8125, 1.2) | 2 = 2.36 x 10~ 3 , and \ip(0, -22.1875, 2.4)| 2 = 1.25 x 10~ 3 . The repulsive force from the 
condensate remaining in the trap, combined with free expansion under gravity, creates a 'crescent-shaped' profile 
in the output. 



The mutual repulsion between parent and daughter leads to a small recoil of the more massive, trapped 
condensate. The subsequent motion of the parent can be inferred by studying its centre of mass, given by 
rj c = J J i]\ip\ 2 d£dr). Simulations of a condensate in free-fall show that calculations of the centre of mass 
agree with the classical result, rj c = —Gt 2 , to five significant figures. The centre of mass coordinate of 
the parent condensate along with its second time derivative (acceleration), are displayed in Fig. 4. The 
repulsive force from the output pulse reaches a maximum at r ~ 0.75. Subsequently, as the coupling term 
falls to zero due to the decreasing overlap, the parent undergoes harmonic motion at the trap frequency, as 
expected. Fig. 4 also illustrates that the recoil is less for smaller output fractions. Given that the output 
coupling for small / does not seriously perturb the parent, further coherent pulses may be produced, and 
one would expect that the recoil on the parent would be negligible for slow, continuous output coupling. 
Given the agreement between these dynamical simulations and experiment, we surmise that the same nu- 
merical methods can be employed to model the dynamics of the condensate in other interesting situations. 
An example is discussed in the following section. 
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FIG. 4. The vertical position (solid line) and acceleration (dashed) of the centre of mass of the trapped parent 
condensate (C = 200 and G = 3.0) as a function of time. Repulsion from the output pulse induces a recoil, which 
excites the parent into simple harmonic motion. Results are presented for / = 0.20 (thick lines) and / = 0.05 
(thin lines). 



IV. FLOW THROUGH A CONSTRICTION 



We now consider the flow of a weakly-interacting condensate through a constriction. In the 2D simulations 
presented here, the condensate is released from a circular trap and falls under gravity down a focused, 
far-off resonance laser beam. The beam is approximated by a harmonic potential in which the spring 
constant varies with position along z. Any variation in potential along the z-axis (e.g. due to an increase 
in intensity near to the focus) is neglected. This situation may be realized using a blue-detuned doughnut 
mode laser (see e.g. [^9|). The model also assumes that there is no delay between release from the trap 
and switching on the laser beam. Under this model, for r > 0, the evolution of the condensate is given 
by the equation: 



id T ip 



1 + (vl/Vo) 



!+(»?- Vrn) 2 /Vo 



(ii) 



so that the horizontal oscillation frequency, u> — (rfc +Tjm)/[ 1 lo + iv^Vm) 2 ], is matched to the initial trap 
frequency at the 'release' position, r/ = 0, and increases by a factor of [1 + (j] m /r]o) 2 ] at r\ — r\ m . The 
limit t]q — > oo corresponds to a collimated laser beam. In this case, the simulations show that the laser 
acts as a waveguide for output-coupled condensates. 

An interesting situation arises when the laser is tightly-focused. Fig. 5 displays results of the simulation 
for a small nonlinear coefficient (C = 10). The density plot at r = 2.0 illustrates that three peaks 
are formed above the constriction. Scaling the ID wavefunctions in Section 2 illustrates that the peaks 
approximately match the positions of the maxima in Fig. 2, suggesting that the second excited state is 
populated and therefore parity is conserved. The probability that the condensate will be excited from 
the ground state is related to the rate of change in the Hamiltonian compared to the excitation spectrum. 
For example, in the non-interacting limit (C — > 0) , we can derive an approximate criterion for excitation 



to occur based on adiabaticity: 4G 1 / 2 



Tjn 



11/2 



/ 770 3> A/i, where A/i is the difference in energy between 



the ground and excited states. Thus, a small constriction creates a sudden change in the Hamiltonian 
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and mode excitation occurs. Conversely, for wider constrictions, the Hamiltonian changes slowly and the 
fluid can adjust adiabatically. 
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FIG. 5. Flow of a weakly-interacting condensate (C = 10, G = 3.0) through a constriction, where the oscillation 
frequency increases by a factor of 11.5 between rj = 0.0 and rj m = —6.0. Contours of the density are shown at 
r = 2.0 (offset left) and r = 2.5 (offset right). The excited 'tail' of the condensate remains above the constriction 
while the ground state falls under gravity, leading to spatial separation of the eigenstates. The contours are drawn 
with an equal spacing of 7 x 10 -3 . 



At subsequent times (r = 2.5 in Fig. 5) the excited 'tail' remains above the constriction, while the 'bulk' of 
the condensate in the ground state continues to fall. To understand this separation of eigenmodes, consider 
the simple case of a non-interacting condensate, C = 0. As the condensate falls into the constriction, the 
increase in oscillation frequency, Slo, leads to a rise in the chemical potential by 5fi = (n + 1/2) Slo. This 
increase in energy is provided by a change in the gravitational potential, Grj. Equating these energies 
gives an expression for the equilibrium position, r\ e : 




Vn 



G 

So, for the n-th mode to be trapped, the following inequality must hold: 



n>2G[r] 



'If, 



Vo 



1 



(12) 



(13) 



In the case of Fig. 5, equation ( |l2] ) predicts rj e ~ —4.3 for n — 2, just below the actual position as expected, 
because the small nonlinearity tends to push the mode upwards. In addition, the inequality ( |l3| ) gives 
n > 1.2, confirming that only excited states are trapped above the constriction. The constriction thus acts 
as a mode 'filter', and could provide a convenient means to generate and study the decay of excitations. 
However, for a light-induced potential this would require tight-focusing: typically spot sizes of 1 /im or 
less. 

For large C, the interaction term in the chemical potential dominates. Thus, in addition to the rise in 
energy due to the changing uj, there is an increase in the density and therefore the interaction energy. As 
a result, the relative difference between eigenvalues is much lower than for C = 0, and mode separation 
tends not to occur. In Fig. 6, one can see that, for C — 200, two maxima form above a wide constriction. 
Parity is still conserved, but the 'tail' in this case consists mainly of a superposition of ground and second 
excited states. If narrower constrictions are used, many modes are excited, and the situation becomes 
increasingly complex. 
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FIG. 6. Flow of a condensate (C = 200, G = 3.0) through a constriction at a time r = 2.0, where the oscillation 
frequency increases by a factor of 2.5 between r\ — 0.0 and 77 = —6.0. As for the weakly-interacting case, the 
potential induces excitations in the condensate. Contours are drawn with an equal spacing of 3.03 x 10~ 3 . 

V. SUMMARY 

In this paper, we present ID and 2D solutions of the Gross-Pitaevskii equation. ID wavefunctions for 
the first and second excited states are presented for various interaction strengths, and provide useful 
information for the subsequent interpretation of our 2D simulations. 2D ground state solutions are also 
calculated for a circular geometry using a finite difference method. We note that convergence towards 
the interaction dominated Thomas-Fermi limit is slower in two than in one dimension. 
The 2D wavefunctions were used as initial conditions in the following 2D simulations. Firstly, an output 
coupler for magnetically-trapped Bose condensed atoms was simulated by solution of a pair of coupled 
Gross-Pitaevskii equations, to represent the output pulse and parent condensate. We demonstrate that a 
combination of atomic interactions and gravity leads to curvature in the output, in qualitative agreement 
with experiment. The simulations also demonstrate that dipolar centre-of-mass modes are excited in 
the remaining trapped condensate. The amplitude of the mode is particularly sensitive to the output 
fraction, and could provide a useful diagnostic in experiments. Secondly, we presented simulations of 
the flow of a condensate dropped under gravity into a tightly-focused laser beam. We demonstrate that 
a tight constriction may be used to excite and separate eigenmodes in a weakly-interacting condensate. 
This could provide a useful technique to study the decay of excited condensate modes. 
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